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1.  Introduction 


Prediction  of  the  in-flight  motion  of  projectiles  requires  the  determination  of  the  aerodynamic 
forces  and  moments  that  act  on  the  body  during  the  course  of  its  trajectory.  These  aerodynamic 
forces  and  moments  may  be  determined  from  free-flight  aerodynamics  ranges  (7),  wind  tunnel 
experiments  (2),  or  theoretical  approaches,  such  as  computational  fluid  dynamics  (CFD).  In  the 
supersonic  flight  regime,  static  aerodynamic  coefficients  such  as  drag  and  pitching  moment  have 
been  routinely  computed  using  CFD  for  a  number  of  years  (3, 4).  Much  less  routine  are 
predictions  of  dynamic  aerodynamic  coefficients,  such  as  the  pitch  damping  and  Magnus 
coefficients.  A  computational  capability  for  predicting  the  Magnus  forces  and  moments  for 
smooth  medium-  and  large-caliber  flight  bodies  has  existed  since  the  early  1980s  (3, 4).  With 
die  recent  development  of  capabilities  for  predicting  the  damping  derivatives  (J),  it  is  now 
possible  to  predict  the  stability  characteristics  and  free-flight  motion  for  projectiles  using  data 
that  are  derived  solely  from  CFD. 

While  previous  research  work  has  focused  on  predictions  of  the  individual  aerodynamic 
coefficients,  this  report  focuses  on  the  predictions  of  aerodynamic  performance,  in  particular, 
aerodynamic  stability.  This  report  attempts  to  demonstrate,  through  a  particular  case  study,  the 
suitability  and  accuracy  of  Navier-Stokes  computational  techniques  for  the  prediction  of 
projectile  performance  and  stability.  Until  recently,  this  was  not  possible  since  methods  had  not 
been  implemented  to  predict  all  the  necessary  coefficients. 

As  a  demonstration  of  the  capability,  this  report  will  present  results  for  a  family  of  axisymmetric 
projectiles  in  supersonic  flight.  The  particular  configuration  selected  for  this  computational 
study  has  been  extensively  tested  in  aeroballistic  ranges,  and  high-quality  experimental  data  have 
been  obtained  (7).  The  configuration  consists  of  5-,  7-,  and  9-caliber  length  versions.  By 
varying  the  mass  distribution,  three  different  centers  of  gravity  (CG)  were  tested  for  each  body 
length.  A  schematic  of  the  configuration  is  shown  in  Figure  1. 

The  computational  results  were  obtained  by  applying  a  parabolized  Navier-Stokes  (PNS) 
approach  to  compute  the  attached  viscous  flow  over  the  forebody  of  the  projectile.  Aerodynamic 
forces  and  moments  were  then  extracted  from  the  computed  flow  fields  by  integration  of  the 
predicted  pressure  and  shear  stresses  acting  on  the  projectile  surface.  To  compute  the 
recirculating  flow  at  the  base  of  the  projectile  and  subsequently  the  base  drag,  a  time-marching 
Navier-Stokes  approach  was  applied. 

Using  the  predicted  aerodynamics  coefficients,  parameters  that  characterize  the  in-flight  motion 
are  subsequently  evaluated,  including  the  gyroscopic  and  dynamic  stability  factors  and  the 
projectile’s  fast  and  slow  mode  frequencies  and  damping  coefficients.  These  parameters  are  then 
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Figure  1 .  Schematic  of  Army-Navy  Spinner  rocket  (ANSR)  configuration. 


used  to  predict  the  free-flight  motion  of  the  projectile.  In  each  case,  the  computational  approach 
is  validated  by  comparison  with  experimental  data. 


2.  Computational  Approach 


The  supersonic  viscous  flow  field  about  the  projectile  configuration  was  obtained  by  applying  a 
parabolized  Navier-Stokes  technique  to  determine  the  attached  flow  field  over  the  projectile 
forebody  and  a  time-marching  Navier-Stokes  technique  to  determine  the  recirculating  flow  in  the 
base  region.  Both  of  these  numerical  techniques  are  briefly  described  in  the  following 
subsections. 

2.1  The  PNS  Approach 

Computation  of  the  attached  viscous  flow  field  about  the  projectile  forebody  was  accomplished 
by  solving  the  thin-layer  Navier-Stokes  equations  using  the  PNS  technique  of  Schiffand  Steger 
(<5).  Using  the  PNS  technique,  computational  results  are  obtained  by  marching  through  the  grid 
from  the  projectile  nose  to  the  base.  This  technique  is  applicable  in  the  supersonic  flow  regime 
and  requires  that  the  flow  field  contain  no  regions  of  flow  separation  in  the  axial  direction. 
Because  the  computational  approach  requires  only  a  single  sweep  through  the  computational 
grid,  it  is  very  efficient,  compared  with  time-marching  approaches  that  require  many  sweeps 
through  the  grid. 

As  is  standard  practice,  prediction  of  the  static  aerodynamic  coefficients,  such  as  pitching 
moment  and  normal  force  coefficients,  are  performed  by  computing  the  flow  field  about  the 
projectile  at  a  constant  angle  of  attack.  The  Magnus  coefficients  are  obtained  from  a 
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computation  of  the  flow  field  at  a  constant  angle  of  attack  and  spin  rate.  Since  the  static 
aerodynamics  for  spinning  axisymmetric  bodies  are  generally  invariant  with  spin  rate,  it  is 
possible  to  obtain  the  static  aerodynamic  coefficients  and  the  Magnus  coefficients  from  a  single 
computation.  The  pitch-damping  coefficients  are  determined  from  the  side  force  and  moment 
acting  on  the  projectile  undergoing  steady  coning  motion.  Steady  coning  motion  is  defined  as 
the  motion  performed  by  a  missile  flying  at  a  constant  angle  with  respect  to  the  freestream 
velocity  vector  and  undergoing  a  rotation  at  a  constant  angular  velocity  about  a  line  parallel  to 
the  freestream  velocity  vector  and  coincident  with  the  projectile’s  CG.  Coning  motion  is,  in  fact, 
a  specific  combination  of  two  orthogonal  planar  pitching  motions  plus  a  spinning  motion.  From 
linear  flight  mechanics  theory,  it  can  be  shown  that  the  pitch-damping  coefficients  are  related  to 
the  side  force  and  moment  due  to  steady  coning  motion  (5). 

The  flow  field  predictions  of  the  projectile  in  steady  coning  motion  have  been  performed  in  a 
novel  rotating  coordinate  frame  that  rotates  at  the  roll  rate  or  coning  rate  of  the  projectile.  The 
fluid  flow  relative  to  the  rotating  coordinate  frame  does  not  vary  with  time,  allowing  the  steady 
(non-time  varying)  Navier-Stokes  equations  to  be  applied.  To  implement  the  rotating  coordinate 
frame,  the  governing  equations  have  been  modified  to  include  the  effect  of  centrifugal  and 
Coriolis  forces.  The  steady,  thin-layer  Navier-Stokes  equations  are  shown  in  equation  1: 

dE  d¥  dG  «  1  as 

d£  dT]  d£  Re  dg 

Here,  E ,  F ,  and  G  are  the  inviscid  flux  vectors,  S  is  the  viscous  flux  vector,  and  H  is  the 
source  term  containing  the  Coriolis  and  centrifugal  force  terms,  which  result  from  the  rotating 
coordinate  frame.  Each  of  these  matrices  is  a  function  of  the  dependent  variables  represented  by 
the  vector  q (p,  pu,  p\,  />w,  e),  where  p  and  e  are  the  density  and  the  total  energy  per  unit  volume, 
and  u,  v,  and  w,  are  the  velocity  components  in  x,  y,  and  z  directions.  Further  details  of  the  flux 
terms  and  the  source  term  containing  the  Coriolis  and  centrifugal  forces  are  given  by  Schiff  and 
Steger  (6)  and  Weinacht  et  al.  (5),  respectively. 

The  thin-layer  equations  are  solved  using  the  PNS  technique  of  Schiff  and  Steger  (6).  Following 
the  approach  of  Schiff  and  Steger,  the  governing  equations,  which  have  been  modified  here  to 
include  the  Coriolis  and  centrifugal  force  terms,  are  solved  using  a  conservative,  approximately 
factored,  implicit  finite-difference  numerical  algorithm  as  formulated  by  Beam  and  Wanning  (7). 
A  fully  turbulent  boundary  layer  has  been  modeled  using  the  Baldwin-Lomax  turbulence  model 
(8). 

The  computations  presented  here  were  performed  using  a  shock-fitting  procedure  reported  by 
Rai  and  Chaussee  (9).  This  procedure  solves  the  five  Rankine-Hugoniot  jump  conditions,  two 
geometric  shock-propagation  conditions,  and  one  compatibility  equation  to  determine  the  values 
of  the  five  dependent  variables  immediately  behind  the  shock,  as  well  as  the  position  of  the 
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shock.  The  shock-fitting  procedure  was  appropriately  modified  when  coning  motion  was 
employed  (5). 

The  computational  results  presented  here  were  obtained  using  a  grid  that  consisted  of  60  points 
between  the  body  and  the  shock.  In  the  circumferential  direction,  gridding  was  performed  over 
the  entire  body  360°  using  36  points  spaced  at  10°  increments.  Approximately  75  axial 
marching  steps  were  used  for  each  caliber  of  body  length.  Because  several  of  the  aerodynamic 
coefficients  are  strongly  influenced  by  viscous  effects,  the  boundary  layer  was  adapted  by 
monitoring  the  nondimensional  boundary  layer  coordinate,  y+,  along  the  leeside  of  the  body  and 
adjusting  the  normal  grid  spacing  so  that  the  y+  was  ~3.  This  approach  has  been  successfully 
applied  in  previous  studies  (3). 

The  computations  have  been  performed  on  a  variety  of  computer  platforms.  Typical  solutions 
for  a  single  set  of  conditions  (Mach  number,  coning  rate,  angle  of  attack,  etc.)  required  5-10  min 
total  CPU  time  on  a  Cray  SV1  computer. 

The  selection  of  grid  parameters  for  the  current  study  is  consistent  with  those  used  in  the 
previous  studies  that  focused  on  the  individual  aerodynamics  coefficients  (5-5).  The  grid 
requirements  shown  in  these  previous  studies  indicate  that  accurate  results  for  the  complete  set  of 
aerodynamics  coefficient  can  be  obtained  with  a  single  grid.  As  part  of  the  current  effort,  a  grid 
refinement  study  for  the  length-to-diameter-ratio  ([L/D]  =  9)  body,  middle  CG  position  was  also 
performed.  A  more  refined  grid  with  90  points  between  the  body  and  the  shock,  72 
circumferential  points  and  150  axial  marching  steps  for  each  caliber  of  body  length  was  utilized. 
A  comparison  of  the  results  obtained  with  the  baseline  and  refined  grids  indicated  a  maximum 
difference  of  <1.5%  between  the  two  complete  sets  of  aerodynamic  coefficients. 

2.2  The  Time-Marching  Approach 

For  the  recirculating  flow  in  the  base  region  of  the  flow  field,  the  PNS  technique  can  no  longer 
be  used,  and  a  more  generally  applicable  approach  such  as  the  time-marching  technique  must  be 
applied.  The  time-marching  technique  is  based  on  the  unsteady  or  time-dependent  Navier- 
Stokes  equations.  For  steady  flow  applications,  such  as  those  discussed  here,  the  solution  is 
obtained  by  marching  an  initial  “guessed”  solution  to  a  final  converged  solution  in  a  time- 
iterative  manner.  This  requires  the  entire  flow  field  to  be  updated  at  each  time  step  and 
numerous  time  steps  are  required  as  the  solution  evolves  to  the  converged  steady-state  result. 
Because  of  this,  the  time-marching  approach  is  much  more  computationally  intensive  in  terms  of 
its  run-times  and  computer  memory  requirements  than  the  PNS  technique. 

For  the  predictions  presented  here,  the  time-dependent  approach  was  applied  to  compute  the 
axisymmetric  (zero  angle  of  attack)  base  flow  as  a  means  of  evaluating  the  zero-yaw  drag.  The 
small,  but  measurable,  effect  of  yaw  on  drag  was  not  considered  in  this  study.  The  effect  of  the 
base  flow  on  the  other  coefficients  was  also  assumed  to  be  small. 
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The  particular  time-marching  technique  applied  here  is  the  implicit,  partially  flux-split,  upwind 
numerical  scheme  developed  by  Steger  and  his  colleagues  (10, 11).  The  technique  is  based  on 
the  flux-splitting  approach  of  Steger  and  Warming  (12)  and  is  often  referred  to  as  the  F3D 
technique.  The  upwind  scheme  provides  natural  numerical  dissipation,  better  stability  properties 
and  better  computational  efficiency  compared  with  implicit  central  differencing  schemes 
(10-12). 

The  technique  is  based  on  the  time-dependent  thin-layer  Navier-Stokes  equations  shown  in 
equation  2: 

9q  9E  9F  5G  «  1  dS  ... 

—+ —  +  —  +  —  +  H  = - ;  (2) 

at  a£  drj  d£  Re  d£ 

The  flux-splitting  algorithm  of  the  F3D  approach  reduces  the  computational  requirements  by 
eliminating  one  of  the  three  implicit  inversions  associated  with  the  central  difference  Beam- 
Warming  algorithm  (7).  The  two-factor  implicit  algorithm  involves  two  sweeps  through  the  grid 
at  each  time  step.  The  first  sweep  involves  inverting  the  block  tridiagonal  system  of  equations  in 
the  ^-direction  along  grid  lines  of  constant  £  and  tj  to  determine  the  intermediate  solution 
variable.  During  the  second  sweep,  a  second  block  tridiagonal  system  of  equations  is  inverted  in 
the  7-direction  along  grid  lines  of  constant  £and  ^to  determine  the  dependent  variable. 

The  particular  version  of  the  code  applied  here  has  a  multizone  capability  that  allows  the 
computational  grid  to  be  decomposed  into  smaller  blocks.  This  allows  greater  flexibility  in 
gridding  of  flight  bodies,  particularly  when  it  is  desirable  to  preserve  sharp  comers.  The 
multizone  capability  also  allows  the  computer  memory  requirement  for  the  computation  to  be 
reduced  since  only  one  zone  needs  to  be  in  memory  at  a  time  and  the  other  zones  can  be  stored 
on  disk.  This  capability  was  implemented  by  Sahu  and  Steger  (11). 

The  zonal  gridding  approach  was  used  to  compute  the  recirculating  viscous  flow  in  the  base 
region.  A  large  grid  zone  was  located  in  the  base  region  and  extended  axially  from  the  projectile 
base  6.35  calibers  down  stream.  In  the  radial  direction,  the  grid  extended  outward  from  the  axis 
of  symmetry  to  a  radial  location  3.5  calibers  above  the  base  comer. 

To  allow  the  flow  to  properly  expand  around  the  base  comer,  a  small  zone,  one  body  diameter  in 
length,  was  placed  on  the  forebody  upstream  from  the  base  comer.  This  grid  had  a  one-plane 
overlap  with  the  base  grid  and  a  one-to-one  mapping  of  the  grid  points  in  the  overlap  region. 

The  dimensions  of  the  base  grid  were  120  axial  points,  50  radial  points  along  the  base,  and  50 
radial  points  above  the  base  comer.  The  forebody  grid  had  30  axial  points  and  50  radial  points. 
Clustering  of  the  grid  points  in  the  radial  direction  was  performed  to  resolve  the  boundary  layer 
on  the  forebody  and  the  shear  layer  in  the  base  region.  The  algebraic  turbulence  model  of  Chow 
(13)  was  applied  in  the  base  region. 
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3.  Results 


As  a  first  step  in  evaluating  the  aerodynamic  performance  of  a  flight  vehicle,  the  aerodynamic 
force  and  moment  coefficients  must  be  evaluated.  The  aerodynamic  coefficients  allow  forces 
and  moments  due  to  complicated  three-dimensional  (3-D)  motions  to  be  characterized  as  a 
combination  of  aerodynamic  forces  and  moments  due  to  simple  specific  motions.  The  individual 
aerodynamic  coefficients  can  provide  insight  into  the  aerodynamic  performance  of  the  flight 
vehicles. 

3.1  Aerodynamic  Moments 

For  symmetric  missiles,  the  aerodynamic  moments  are  typically  modeled  using  the  following 
moment  expansion  which  is  cast  in  a  non-rolling  coordinate  frame  (14).  The  moment 
formulation  makes  use  of  complex  variables  to  separate  the  moment  components,  Cm  and  Cn 

which  produce  rotations  in  the  vertical  and  horizontal  planes,  respectively.  The  moment  about 
the  longitudinal  axis  of  the  vehicle,  Ci,  can  be  treated  separately. 


Cm+iC„  = 

\V  ) 

(3) 

c.-feev 

IvJ  'p 

(4) 

In  the  moment  formulation,  the  pitching  moment  coefficient,  Cm^ ,  produces  a  moment 

proportional  to  the  complex  angle  of  attack,  g  .  For  small  angle-of-attack  flight,  the  cosine  of 
the  total  angle  of  attack,  y ,  is  nearly  1.  This  allows  the  pitch-damping  moment  coefficient, 

c"v,+c"»  ,  which  produces  moments  proportional  to  the  angular  rate  <*' ,  to  be  treated  as  a 
single  coefficient.  The  Magnus  moment  coefficient,  C„  ,  accounts  for  a  side  moment  due  to 

“par 

flow  asymmetries  from  a  combination  of  spin  and  angle  of  attack. 

Over  the  past  decade  and  a  half,  CFD  techniques  have  been  developed  to  predict  the  various 
aerodynamic  coefficients  previously  discussed.  The  pitching  moment  coefficient,  Cma ,  has 

been  routinely  determined  for  various  configurations  and  flight  regimes  for  well  over  a  decade. 
The  prediction  of  the  Magnus  coefficients  for  axisymmetric  projectiles  at  supersonic  velocities 
was  demonstrated  in  the  early  1980s  (3).  The  Magnus  problem  is  still  a  research  topic  for 
complicated  vehicles  and  at  transonic  velocities.  More  recently,  techniques  for  predicting  the 
pitch-damping  coefficient,  C,^  +  Cm. ,  based  on  the  solution  of  the  Navier-Stokes  equations 

have  been  developed  (5).  Although  the  development  of  a  capability  to  predict  each  individual 
aerodynamic  coefficient  can  be  considered  an  important  research  accomplishment,  it  was  not 
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possible  to  predict  the  stability  characteristics  of  projectiles  using  CFD  until  a  complete 
capability  for  predicting  all  the  aerodynamic  coefficients  was  developed.  The  prediction  of  the 
pitch  damping  coefficient  represented  that  final  step. 

Figures  2-4  show  results  of  CFD  predictions  for  the  pitching  moment,  Magnus  moment  and 
pitch-damping  moment  coefficients  performed  as  part  of  the  current  study.  Results  for  a 
ffee-flight  Mach  number  of  1 .8  are  shown  for  each  of  the  three  CG  locations  and  body  lengths. 

In  general,  the  CFD  predictions  appear  to  be  in  good  agreement  with  the  experimental  data 
obtained  from  range  firings  of  these  projectiles.  Results  obtained  at  Mach  2.5  also  show  similar 
agreement  with  the  range  data. 

The  aerodynamic  moment  about  the  longitudinal  axis  of  an  axisymmetric  vehicle  can  be 
modeled  using  a  single  aerodynamic  coefficient.  The  roll-damping  coefficient,  C]p ,  represents 

the  aerodynamic  moment  which  opposes  the  spinning  motion  of  the  vehicle.  For  axisymmetric 
flight  bodies,  the  roll-damping  coefficient  is  produced  purely  by  the  viscous  stresses  acting 
circumferentially  on  the  body  surface.  From  a  modeling  standpoint,  proper  prediction  of  the  roll¬ 
damping  coefficient  requires  adequate  grid  resolution  of  the  boundary  layer.  A  comparison  of 
the  predicted  roll-damping  coefficient  with  range  data  is  shown  in  Figure  5.  Both  the 
computational  predictions  and  the  range  data  show  a  nearly  linear  increase  in  the  roll-damping 
coefficient  with  increasing  body  length.  The  predicted  result  shows  an  overprediction  of  ~10%. 
Similar  results  were  observed  at  Mach  2.5.  The  predictions  showed  a  <1%  variation  in  the  roll¬ 
damping  coefficient  at  angles  of  attack  to  6°. 

3.2  Aerodynamic  Forces 

The  aerodynamic  forces  can  also  be  described  in  a  manner  that  is  analogous  to  the  aerodynamic 
moments.  The  transverse  force  components  CY  and  Cz  act  along  the  y  and  z  axes  of  the  non¬ 
rolling  coordinate  frame  and  consist  of  components  proportional  to  the  complex  yaw  £  and 
yawing  rate  £ '  •  The  third  force  component,  Cx ,  which  acts  along  the  axial  direction,  is  loosely 
coupled  to  the  transverse  forces  and  can  be  treated  separately.  This  force  component  consists 
predominantly  of  the  drag  force. 


The  normal  force  coefficient  CNa  acts  normal  to  the  longitudinal  axis  of  the  body  in  the  pitch- 
plane.  The  force  expansion  also  contains  a  damping  force  coefficient,  CN  +  CN .  ,  that  produces 

Q  or 

a  force  proportional  to  the  yawing  rate  £ ' .  The  Magnus  force  coefficient  CYpa  accounts  for  the 
side  force  produced  by  a  combination  of  spin  and  yaw. 
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Figure  2.  Pitching  moment  coefficient  vs.  CG  position,  Mach  =  1.8,  ANSR. 


Figure  3.  Magnus  moment  coefficient  vs.  CG  position,  Mach  =  1.8,  ANSR. 
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Figure  4.  Pitch-damping  moment  coefficient  vs.  CG  position,  Mach  =  1.8,  ANSR. 
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Figures  6  and  7  show  comparisons  of  the  predicted  and  experimentally  determined  normal  force 
and  Magnus  force  coefficients  at  Mach  1 .8  as  a  function  of  body  length.  The  predicted  normal 
force  coefficient  is  predicted  to  within  the  scatter  of  the  range  data  and  shows  little  influence  of 
body  length.  Both  the  predicted  and  experimentally  determined  Magnus  force  coefficient 
increase  with  body  length.  The  Magnus  force  determined  from  the  variation  of  the  Magnus 
moment  with  CG  location  is  -30%  less  than  the  predictions  for  the  L/D  =  7  and  L/D  =  9  bodies. 
The  Magnus  force  coefficients  determined  from  the  swerving  (CG)  motion  for  individual  shots 
shows  significant  scatter,  indicating  the  difficulty  in  accurately  determining  this  coefficient.  The 
pitch-damping  force  coefficient  is  not  shown  as  it  does  not  influence  the  projectile  motion  in  a 
measurable  manner.  Comparisons  of  predictions  with  experimental  data  are  presented  by 
Weinacht  et  al.  (J). 
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Figure  6.  Normal  force  coefficient  vs.  body  length,  Mach  =  1.8,  ANSR. 


The  drag  coefficient  is  also  very  important  for  predicting  projectile  flight  performance  as  it  is 
one  of  the  primary  aerodynamic  coefficients  affecting  range.  The  drag  coefficient  is  also 
required  for  the  stability  assessment,  although  its  effect  is  small.  The  total  drag  on  the  projectile 
is  composed  of  contributions  from  the  forebody  and  the  base.  The  forebody  drag  consists  of 
both  pressure  and  viscous  drag  components.  The  base  drag  results  from  the  pressure  acting  on 
the  projectile  base.  The  forebody  drag  can  be  predicted  using  the  PNS  technique.  The 
recirculating  flow  in  the  base  region  is  predicted  using  the  time-marching  technique  from  which 
the  base  drag  can  be  determined. 
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Figure  7.  Magnus  force  coefficient  vs.  body  length,  Mach  =  1.8,  ANSR. 

A  comparison  of  the  predicted  zero-yaw  drag  with  range  results  is  shown  in  Figure  8.  Drag 
predictions  at  Mach  1.8  and  2.5  are  shown  as  a  function  of  body  length.  Both  the  experimental 
results  and  the  computational  results  show  a  slight  increase  in  zero-yaw  drag  with  increasing 
body  length  and  decreasing  Mach  number.  The  computed  drag  coefficient  is  predicted  to  within 
the  accuracy  of  the  range  data.  The  increase  in  the  drag  with  body  length  is  due  primarily  to  the 
increase  in  the  wetted  surface  of  the  longer  body.  The  predicted  results  showed  only  minor 
variations  in  the  base  drag  with  body  length. 

At  Mach  1.8  for  the  L/D  =  9  body,  the  pressure  and  viscous  forebody  drag  account  for  39%  and 
25%  of  the  total  drag  with  the  base  drag  accounting  for  the  remaining  36%  of  the  drag. 
Predictions  of  the  forebody  pressure  drag  are  within  1%  of  previous  predictions  using  empirical 
(15)  or  inviscid  ( 1 )  approaches. 

An  estimate  of  the  predicted  skin  friction  drag,  CDy ,  can  be  obtained  using  the  relation 
developed  by  Charters  and  Kent  (16). 

CDv  =-4C,p.  (7) 

The  estimated  skin  friction  drag  determined  from  the  predicted  roll-damping  coefficient 
underestimates  the  predicted  skin  friction  drag  by  3%-10%  for  the  predictions  presented  here. 
(This  comparison  gives  an  indication  of  the  accuracy  of  the  Charters  and  Kent  relation.)  The 
estimated  skin  friction  obtained  from  the  range  measurements  of  the  roll-damping  coefficient  is 
13%-22%  less  than  the  predicted  skin  friction. 
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Figure  8.  Zero-yaw  drag  coefficient  vs.  CG  position,  Mach  =  1.8,  ANSR. 


Although  these  estimates  indicate  a  potential  10%  overprediction  of  the  skin  friction  drag,  the 
contribution  to  the  total  drag  would  result  in  an  overprediction  of  only  a  few  percent. 

Design  codes  can  provide  a  more  cost-effective  alternative  for  predicting  the  base  drag,  thus 
eliminating  the  need  for  a  time-dependent  base  drag  computation.  The  semi-empirical  approach 
of  McCoy  (15)  appears  to  overpredict  the  base  drag  by  10%-20%  compared  to  the  CFD  result. 
When  the  estimated  base  drag  is  added  to  the  forebody  drag  CFD  predictions,  the  total  drag  is 
overpredicted  by  <9%  compared  to  the  CFD  results.  For  stability  computations,  where  the  drag 
coefficient  plays  a  minor  role,  acceptable  accuracy  can  be  obtained  using  base  drag  estimates. 

3.3  Stability  and  Free-Flight  Motion 

Once  the  aerodynamic  forces  and  moments  have  been  determined,  the  equations  of  motion  can 
be  integrated  and  the  projectile’s  motion  can  be  determined  for  various  initial  conditions. 
Because  the  equations  of  motion  are  nonlinear,  numerical  integration  is  usually  required  to 
predict  the  projectile  motion  from  launch  to  impact.  These  nonlinearities  are  typically  weak,  and 
for  portions  of  the  flight  trajectory  the  equations  can  be  integrated  and  closed-form  solutions  for 
the  projectile  motion  can  be  obtained. 

For  a  symmetric  missile,  the  closed-form  solution  for  the  angular  response  as  a  function  of  the 
distance  downrange,  s,  takes  the  following  form  (14): 
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£  =£g  +Kj  +K2  «  fi  +  ia 

Kx  =  K10  exp[i(^10  +  A's)  +  ^s]  (8) 

K2  =  K20  exp[i(^20  +  f2  s)+  ^s]. 

The  complex  variable,  £  ,  consists  of  the  vertical  and  horizontal  components  of  the  angle  of 
attack,  a  and  /?.  The  angular  motion  contains  a  gravity  term,  £g ,  which  is  small  and,  for  the 

remainder  of  the  discussion,  this  term  will  be  ignored.  The  second  and  third  terms  of  equation  8 
are  expressions  of  the  amplitudes  of  two  components  of  the  angular  motion,  referred  to  as  the 
fast  and  slow  modes,  respectively.  The  complex  angle  of  attack,  £  ,  can  be  displayed 
graphically  along  with  the  fast  and  slow  mode  arms,  Ki  and  K2.  Graphically,  the  solution  is 
shown  by  Figure  9  and  is  referred  to  as  a  damped  epicycle.  It  is  helpful  to  regard  this  figure  as 
the  path  that  the  nose  of  the  projectile  would  follow  as  it  flies  downrange.  The  fast  and  slow 
mode  frequencies  of  the  motion,  </>{  and  <f>'2 ,  and  the  damping  rates,  \  and  ,  are  functions  of 

the  aerodynamic  coefficients,  projectile  spin  rate,  and  inertial  properties  of  the  body.  Generally 
speaking,  the  frequencies  are  influenced  by  the  pitching  moment  coefficient,  while  the  damping 
is  influenced  by  the  Magnus  and  pitch  damping  coefficients.  The  equation  also  contains  the 
initial  amplitudes,  K10  and  K2o,  and  phase  angles,  ^10  and  $20,  which  result  from  applying  the 
initial  conditions  during  the  integration  of  the  equations. 


Figure  9.  Projectile  angular  motion  epicycle. 
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Using  the  predicted  aerodynamic  coefficients,  the  fast  and  slow  mode  damping  and  frequencies 
were  computed  for  the  various  body  lengths,  CG  positions,  launch  Mach  numbers,  and  spin  rates. 
Figure  10  shows  the  predicted  slow  mode  frequency  plotted  against  the  experimentally 
determined  value.  (The  results  are  plotted  in  this  fashion  due  to  the  large  number  of  parameters 
considered  here.)  Points  that  fall  close  to  the  solid  line  indicate  that  the  predicted  value  is  close 
to  the  experimentally  determined  value.  There  is  good  agreement  between  prediction  and 
experiment.  The  slow  mode  frequency  represents  the  critical  comparison  for  the  frequencies 
because  the  fast  mode  frequency  is  dominated  by  the  spin  rate,  which  is  essentially  an  input 
parameter. 


Figure  10.  Comparison  of  computed  slow  mode  frequency  with  measured  slow  mode 
frequency,  ANSR. 


Figures  1 1  and  12  show  comparisons  of  the  predicted  slow  and  fast  mode  damping  of  each  of  the 
bodies.  The  data  are  scattered  about  the  solid  line,  indicating  that  the  prediction  is  generally  in 
good  agreement  with  the  experiment.  No  consistent  bias  between  the  computation  and 
experimental  data  is  seen.  This  indicates  that  there  is  no  consistent  underprediction  or 
overprediction  of  the  damping  rates.  Much  of  the  scatter  in  the  data  is  attributed  to  the 
uncertainty  in  the  experimental  data  rather  than  inaccuracies  in  the  computational  approach. 

Generally,  from  a  design  standpoint,  predicting  and  controlling  damping  rates,  rather  than  the 
frequencies,  is  a  primary  consideration.  One  desirable  requirement  is  that  the  damping  rates,  Jtf 
and  X.2,  be  negative  so  that  the  angular  motion  does  not  grow  during  the  projectile’s  flight. 
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COMPUTED  X 


Figure  1 1 .  Comparison  of  computed  fast  mode  damping  with  measured  fast  mode 
damping,  ANSR. 


Figure  12.  Comparison  of  computed  slow  mode  damping  with  measured  slow  mode 
damping,  ANSR. 
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By  algebraically  manipulating  the  equations  for  the  damping  rates,  three  necessary  conditions 
can  be  derived,  which  will  ensure  that  the  damping  rates  are  negative.  Two  of  these  conditions 
are  shown  in  equation  9: 

sg  >i;  sg  =  f(p,cmc) 

s,<Sd(2-sd);  sd=f(C„t„,CrI,q+Crnn).  (9) 

The  first  condition  states  that  the  gyroscopic  stability  factor,  sg,  must  be  >1  for  the  projectile  to 
be  gyroscopically  stable.  The  second  condition  also  places  an  additional  constraint  on  the 
gyroscopic  stability  factor  through  the  dynamic  stability  factor,  s<j.  As  shown  in  these  equations, 
the  gyroscopic  stability  and  dynamic  stability  factors  are  functions  of  the  aerodynamic 
coefficients.  Additionally,  the  gyroscopic  stability  factor  is  also  a  function  of  the  projectile  spin 
rate,  p.  For  typical  spin-stabilized  projectiles  that  are  statically  unstable,  the  two  stability  criteria 
can  be  satisfied  by  adjusting  the  spin  rate  of  the  projectile,  as  long  as  the  dynamic  stability  factor 
is  between  0  and  2.  If  the  dynamic  stability  factor  is  outside  this  bound,  the  projectile  cannot  be 
stabilized  using  spin.  A  third  condition  for  stability  essentially  requires  that  the  pitch  damping 
coefficient,  Cmq  +  Cmd  to  be  negative.  This  condition  is  rarely  violated  in  practice. 

The  constraints  on  the  gyroscopic  stability  and  dynamic  stability  factors  discussed  previously  can 
be  shown  graphically  on  what  is  often  termed  a  “stability  plot.”  Figures  13  and  14  show  plots  of 
the  results  as  measured  in  the  aerodynamics  range  and  obtained  from  the  computation.  Shown 
on  the  plot  are  results  for  the  each  of  the  three  CG  locations  for  the  three  different  length  bodies, 
various  spin  rates,  and  for  Mach  numbers  of  1 .8  and  2.5.  For  each  firing,  experimental  and 
computational  values  of  the  gyroscopic  stability  factor,  sg,  and  the  dynamic  stability  factor,  s* 
were  computed.  It  is  important  to  realize  that  for  a  given  set  of  conditions  (body  length,  spin  rate, 
launch  Mach  number),  the  computational  results  will  yield  a  unique  value.  However,  because  of 
the  scatter  associated  with  the  experimental  data,  each  firing  yields  a  slightly  different  value. 
Thus,  the  experimental  results  give  the  appearance  of  having  more  data.  It  is  important  to  realize 
this  when  comparing  the  two  sets  of  plots. 

The  computational  and  experimental  results  presented  in  these  figures  show  that  this  family  of 
projectiles  is  statically  unstable  (sg  >  0),  but  has  been  made  gyroscopically  stable  by  appropriate 
selection  of  the  spin  rate  (sg  >1).  The  results  also  show  that  all  of  the  L/D  =  5  projectiles  are 
dynamically  stable.  (The  region  of  dynamic  stability  is  denoted  by  the  shaded  area  on  these 
plots.)  The  results  indicate  that  several  of  the  L/D  =  7  and  L/D  =  9  projectiles  are  marginally 
stable,  though  these  can  be  stabilized  by  selection  of  a  higher  spin  rate  because  the  dynamic 
stability  factor  is  between  0  and  2.  Several  of  the  L/D  =  9  bodies  are  dynamically  unstable  and 
cannot  be  stabilized  by  changing  the  spin  rate.  Further  examination  of  the  data  revealed  that 
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Figure  13.  Stability  plot  for  ANSR  from  range  data. 


Figure  14.  Stability  plot  for  ANSR  from  CFD  predictions. 
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the  IVD  =  9  bodies  which  could  not  be  stabilized  all  had  the  CG  located  in  the  afbnost  position. 
It  is  important  to  note  that  the  computational  approach  correctly  predicts  the  stability 
characteristics  of  the  various  configurations  for  the  conditions  of  dynamic  stability,  marginal 
dynamic  stability,  and  dynamic  instability.  This  demonstrates  the  utility  of  the  computational 
approach  for  the  purposes  of  stability  prediction. 


With  the  prediction  of  the  frequencies  and  damping  rates  of  the  angular  motion,  the  angular 
motion  itself  can  be  predicted  and  compared  with  the  motion  measured  in  the  aerodynamics 
range.  Figures  15  and  16  show  the  angular  motion  of  a  5-caliber-long  projectile  launched  at 
Mach  1.8  and  a  7-caliber  projectile  launched  at  Mach  2.5.  These  3-D  plots  display  the  vertical 
and  horizontal  components  of  angle  of  attack,  a  and  /?,  as  a  function  of  the  distance  downrange, 
Z.  The  thickened  curve  displays  the  3-D  motion.  Two-dimensional  projections  of  the  motion  in 
die  vertical  and  horizontal  planes  are  also  shown.  It  is  useful  to  regard  the  curves  shown  here  as 
the  path  traversed  by  the  nose  of  the  projectile  as  it  flies  downrange.  The  motion  predicted  using 
the  CFD-derived  aerodynamic  coefficients  is  denoted  by  the  solid  lines,  while  the  motion 
obtained  from  the  range  firings  is  shown  by  the  circular  symbols.  The  agreement  between  the 
two  motions  is  excellent.  Comparisons  of  the  in-flight  motion  at  other  flight  conditions  and  for 
the  other  body  geometries  also  showed  good  agreement  with  the  experimental  data. 


Figure  15.  In-flight  angular  motion,  5-caliber  body,  Mach  1.8. 
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Figure  16.  In-flight  angular  motion,  7-caliber  body,  Mach  2.5. 


4.  Conclusion 


CFD  provides  an  alternate  means  for  determining  the  aerodynamic  performance  of  projectiles. 
This  will  allow  projectiles  to  be  designed  with  a  reduced  requirement  for  manufacturing  models 
and  performing  test  firings.  As  a  test  case,  the  aerodynamic  performance  of  a  family  of 
axisymmetric  projectiles  in  supersonic  flight  was  examined.  The  predictions  showed  good 
agreement  with  range  data  for  individual  aerodynamic  coefficients.  Using  these  coefficients,  the 
frequencies  and  damping  rates  of  the  angular  motion  were  computed.  The  results  showed  good 
correlation  with  the  experimental  data  for  various  body  lengths,  center  of  gravity  locations,  Mach 
numbers  and  spin  rates.  The  aerodynamic  coefficients  were  also  used  to  predict  gyroscopic  and 
dynamic  stability  factors.  The  computational  approach  correctly  distinguished  the  conditions  of 
dynamic  instability,  marginal  dynamic  stability,  and  dynamic  stability.  As  a  final  comparison, 
the  in-flight  angular  motion  for  individual  firings  was  compared  with  range  data  and  showed 
good  agreement.  The  results  are  believed  to  properly  represent  the  capability  for  predicting  the 
linear  aerodynamics  for  the  high-speed  low  angle-of-attack  flight  of  aerodynamically  smooth 
axisymmetric  projectiles.  Additional  computational  effort  (and,  in  some  cases,  additional 
research)  may  be  required  to  establish  the  same  capability  at  subsonic  and  transonic  velocities  or 
for  complex  geometries. 
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List  of  Abbreviations  and  Symbols 


Freestream  speed  of  sound 
Drag  coefficient 
Zero-yaw  drag  coefficient 

Viscous  drag  coefficient 


Net  roll  moment  coefficient 


Roll-damping  moment  coefficient, 


afEBl 

l  vj 


Pitching  moment  coefficient, 


\p^  S*,D 


Pitching  moment  coefficient  slope  with  respect  to  angle  of  attack. 


Pitch-damping  moment  coefficient  slope, 


Pitch-damping  moment  coefficient  slope, 


Cm,  +  Cmd  Pitch-damping  moment  coefficient 
C.  Side  moment  coefficient 


Magnus  moment  coefficient, 


Normal  force  coefficient. 


’  fpvx, 
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Normal  force  coefficient  slope  with  respect  to  angle  of  attack, 


Pitch-damping  force  coefficient  slope, 


ac 


N 


a  D 


ac 


N 


da 


Pitch-damping  force  slope, 


+  C> 


Cy»Cz 


c 


*P  a 


D 


e 

E,F,G 

F 

H 

Kj,K2 

^10»^20 

L 

M 

M 

P 

q 

q 

Re 


s 


Pitch-damping  force  coefficient 
Axial  force  coefficient 

Transverse  force  coefficient  in  y  and  z  direction 
Magnus  force  coefficient 

Reference  diameter 

Total  energy  per  unit  volume 

Flux  vectors  in  transformed  coordinates 

Force 

Source  term  in  Navier-Stokes  equations  due  to  rotating  coordinate  frame 
Amplitude  of  fast  and  slow  yaw  arms 
Initial  amplitude  of  fast  and  slow  yaw  arms 

Body  length 
Moment 

1 

Freestream  Mach  number 
Spin  rate,  as  used  in  roll  equations 

Transverse  angular  rate  of  body,  as  used  in  flight  mechanics  equations 
Vector  of  dependent  variables,  as  used  in  Navier-Stokes  equations 
Reynolds  number,  a^p^D/  //„ 

Distance  downrange  normalized  by  reference  length 
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Sd  Dynamic  stability  factor 

sg  Gyroscopic  stability  factor 


Sref 


Reference  area,  Sref  = 


n  D2 

~T~ 


S  Viscous  flux  vector 

t  Time 

u,  v,  w  V elocity  components  in  the  x,  y,  z  directions 

V  Freestream  velocity 

x,  y,  z  Axial,  horizontal,  and  vertical  coordinates 

y+  Nondimensional  boundary  layer  coordinate 

Z  Dimension  downrange  distance 

Greek  Symbols 
a.  Angle  of  attack 

a  Angular  rate  associated  with  angle  of  attack 

P  Yaw  angle 

y  Cosine  of  total  angle  of  attack 

A, ,  ^  Fast  and  slow  mode  damping  rates 

H,  Viscosity 

Z’V’G  Transformed  coordinates  in  the  Navier-Stokes  equations 

£  Complex  angle  of  attack 

£g  Complex  angle-of-attack  component  induced  by  gravity 

p  Density 

^io>^2o  Phase  angle  of  fast  and  slow  mode  arms 

<f>\ ,  <f>i  Fast  and  slow  mode  frequencies 
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Superscripts 


C) 


Rate  of  change  with  respect  to  time 


Rate  of  change  with  respect  to  space 
Referenced  to  non-rolling  coordinate  frame 


Subscripts 

oo  Quantity  evaluated  at  freestream  conditions 
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